## Heterogeneity Analysis ##

pool_state <- pool %>% mutate(
  state = case_when(
    state == 9 ~ "Uttar Pradesh",
    state == 10 ~ "Bihar",
    state == 19 ~ "West Bengal",
    state == 20 ~ "Jharkhand",
    state == 21 ~ "Odisha",
    state == 23 ~ "Madhya Pradesh"
  )) %>% filter(state == "West Bengal")

# Remove observations that contain missing values for electricity supply dimensions
pool_state <- pool_state %>% dplyr::select(elec_respvars, state, village, hh, age, hhsize, bplaay,
                                           scst, sc, st, wave, gender, religion, education, lnexp,
                                           weights) %>% 
  na.omit()

# With Village Fixed Effects ---------------------------------------------------

# Duration Day (Total hours per day)
model.dd.vfe.duration_day.scst = fixest::feols(duration_day ~ scst*wave + age + gender +
                                                 religion + education + hhsize +
                                                 lnexp + bplaay | village, 
                                               data = pool_state, weights = pool_state$weights)

# Duration Night (Hours in the evening)
model.dd.vfe.duration_night.scst = fixest::feols(duration_night ~ scst*wave + age + gender +
                                                   religion + education + hhsize +
                                                   lnexp + bplaay | village, 
                                                 data = pool_state, weights = pool_state$weights)

# Reliability (Outage days per month)
model.dd.vfe.reliability.scst = fixest::feols(reliability ~ scst*wave + age + gender +
                                                religion + education + hhsize +
                                                lnexp + bplaay | village, 
                                              data = pool_state, weights = pool_state$weights)

# Voltage Fluctuation (Voltage fluctuation days per month)
model.dd.vfe.quality_fluc.scst = fixest::feols(quality_fluc ~ scst*wave + age + gender +
                                                 religion + education + hhsize +
                                                 lnexp + bplaay | village, 
                                               data = pool_state, weights = pool_state$weights)

# Low Voltage (Low Voltage days per month)
model.dd.vfe.quality_low.scst = fixest::feols(quality_low ~ scst*wave + age + gender +
                                                religion + education + hhsize +
                                                lnexp + bplaay | village, 
                                              data = pool_state, weights = pool_state$weights)

# With Household Fixed Effects -------------------------------------------------

# Duration Day (Total hours per day)
model.dd.hfe.duration_day.scst = fixest::feols(duration_day ~ scst*wave + hhsize +
                                                 lnexp + bplaay | hh, 
                                               data = pool_state, weights = pool_state$weights)

# Duration Night (Hours in the evening)
model.dd.hfe.duration_night.scst = fixest::feols(duration_night ~ scst*wave + hhsize +
                                                   lnexp + bplaay | hh, 
                                                 data = pool_state, weights = pool_state$weights)

# Reliability (Outage days per month)
model.dd.hfe.reliability.scst = fixest::feols(reliability ~ scst*wave + hhsize +
                                                lnexp + bplaay | hh, 
                                              data = pool_state, weights = pool_state$weights)

# Voltage Fluctuation (Voltage fluctuation days per month)
model.dd.hfe.quality_fluc.scst = fixest::feols(quality_fluc ~ scst*wave + hhsize +
                                                 lnexp + bplaay | hh, 
                                               data = pool_state, weights = pool_state$weights)

# Low Voltage (Low Voltage days per month)
model.dd.hfe.quality_low.scst = fixest::feols(quality_low ~ scst*wave + hhsize +
                                                lnexp + bplaay | hh, 
                                              data = pool_state, weights = pool_state$weights)

## Model and coefficient lists for village and household fixed effects combined ##
models.vfe.dd_elec_scst <- list(model.dd.vfe.duration_day.scst,model.dd.vfe.duration_night.scst,
                                model.dd.vfe.reliability.scst, model.dd.vfe.quality_fluc.scst,
                                model.dd.vfe.quality_low.scst)

models.hfe.dd_elec_scst <- list(model.dd.hfe.duration_day.scst, 
                                model.dd.hfe.duration_night.scst, model.dd.hfe.reliability.scst, 
                                model.dd.hfe.quality_fluc.scst, model.dd.hfe.quality_low.scst)

# Generate latex tables
esttex(models.vfe.dd_elec_scst, 
       se = "cluster", 
       digits = 3,
       dict = c(scst = "SC/ST", wave = "2018", duration_day = "Daily Supply Hours",
                duration_night = "Evening Supply Hours", reliability = "Monthly Outage Days",
                quality_fluc = "Fluctuating Voltage Days", quality_low = "Low Voltage Days",
                village = "Village", hh = "Household"),
       keep = c(":", "2018"),
       fitstat = c("r2"),
       fixef_sizes = T,
       float = F,
       replace = T,
       file = c("./Manuscript/Tables/dd_elec_scst_state_WB_vfe.tex"))

esttex(models.hfe.dd_elec_scst, 
       se = "cluster", 
       digits = 3,
       dict = c(scst = "SC/ST", wave = "2018", duration_day = "Daily Supply Hours",
                duration_night = "Evening Supply Hours", reliability = "Monthly Outage Days",
                quality_fluc = "Fluctuating Voltage Days", quality_low = "Low Voltage Days",
                village = "Village", hh = "Household"),
       keep = c(":", "2018"),
       fitstat = c("r2"),
       fixef_sizes = T,
       float = F,
       replace = T,
       file = c("./Manuscript/Tables/dd_elec_scst_state_WB_hfe.tex"))


# State-wise analysis results collection ---------------------------------------

state_wb_fits <- models.hfe.dd_elec_scst